/*
 * ComplexMath.cs
 * 
 * Copyright (c) 2003-2005, dnAnalytics. All rights reserved.
*/

using System;

namespace dnAnalytics.Math
{
    ///<summary>Provides trigonometric, logarithmic, and other common mathematical functions for complex types</summary>
    public sealed class ComplexMath
    {
        private ComplexMath() { }

        ///<summary>Return the absolute value of a complex type calculated as the Euclidean norm</summary>
        //Based off Numerical Recipes' Cabs.
        public static double Absolute(Complex value)
        {
            if (value.IsNaN())
            {
                return double.NaN;
            }
            else if (value.IsInfinity())
            {
                return double.PositiveInfinity;
            }

            double real = System.Math.Abs(value.Real);
            double imag = System.Math.Abs(value.Imaginary);
            if (value.Real == 0)
            {
                return imag;
            }
            else if (value.Imaginary == 0)
            {
                return real;
            }
            else if (real > imag)
            {
                double temp = imag / real;
                return real * System.Math.Sqrt(1.0 + temp * temp);
            }
            else
            {
                double temp = real / imag;
                return imag * System.Math.Sqrt(1.0 + temp * temp);
            }
        }

        ///<summary>Calculate the complex argument of a complex type.  Also commonly referred to as the phase.</summary>
        public static double Argument(Complex value)
        {
            return System.Math.Atan2(value.Imaginary, value.Real);
        }

        ///<summary>Return the complex conjugate of a complex type</summary>
        public static Complex Conjugate(Complex value)
        {
            return new Complex(value.Real, -value.Imaginary);
        }

        ///<summary>Return the complex cosine of a complex type</summary>
        public static Complex Cos(Complex value)
        {
            if (value.IsNaN())
            {
                return Complex.NaN;
            }
            return Cosh(new Complex(-value.Imaginary, value.Real));
        }

        ///<summary>Return the complex hyperbolic cosine of a complex type</summary>
        public static Complex Cosh(Complex value)
        {
            if (value.IsNaN())
            {
                return Complex.NaN;
            }
            if (value.IsInfinity())
            {
                return Complex.Infinity;
            }
            return new Complex(System.Math.Cos(value.Imaginary) * System.Math.Cosh(value.Real), System.Math.Sin(value.Imaginary) * System.Math.Sinh(value.Real));
        }

        ///<summary>Return the complex exponential of a complex type</summary>
        public static Complex Exp(Complex value)
        {
            if (value.IsNaN())
            {
                return Complex.NaN;
            }

            double exp = System.Math.Exp(value.Real);

            if (value.Imaginary == 0)
            {
                return new Complex(exp);
            }

            return new Complex(exp * System.Math.Cos(value.Imaginary), exp * System.Math.Sin(value.Imaginary));
        }

        ///<summary>Raise 'leftSide' to the power of 'rightSide'</summary>
        ///<param name="leftSide"><b>Complex</b> value as base</param>
        ///<param name="rightSide"><b>Complex</b> value as exponent</param>
        public static Complex Pow(Complex leftSide, Complex rightSide)
        {
            return Exp(rightSide * ComplexMath.Log(leftSide));
        }


        /// <summary>
        /// Return an approximation to Log(1 + x)
        /// </summary>
        /// <param name="x">The real parameter</param>
        /// <returns>The value of Log(1+x)</returns>
        private static double Log1p(double x)
        {
            if (x < -1.0)
            {
                throw new ArgumentOutOfRangeException();
            }
            else if (x > 1.0E+15)
            {
                return System.Math.Log(x);
            }
            else if (System.Math.Abs(x) < 1.0E-16)
            {
                return x;
            }
            else
            {
                double y = 1.0 + x;
                return System.Math.Log(y) - ((y - 1) - x) / y;
            }

        }

        /// <summary>
        /// Return the value of Log(Abs(z)).
        /// </summary>
        /// <param name="value">The complex number.</param>
        /// <returns>The calculated value.</returns>
        private static double LogAbs(Complex value)
        {

            double xabs = System.Math.Abs(value.Real);
            double yabs = System.Math.Abs(value.Imaginary);

            double max, u;
            if (xabs >= yabs)
            {
                max = xabs;
                u = yabs / xabs;
            }
            else
            {
                max = yabs;
                u = xabs / yabs;
            }

            return (System.Math.Log(max) + 0.5 * ComplexMath.Log1p(u * u));
        }



        ///<summary>Return the complex logarithm of a complex type</summary>
        public static Complex Log(Complex value)
        {
            if (value.IsNaN())
            {
                return Complex.NaN;
            }
            return new Complex(ComplexMath.LogAbs(value), System.Math.Atan2(value.Imaginary, value.Real));
            /* double r = ComplexMath.Absolute(value);
             double i = System.Math.Atan2(value.Imaginary, value.Real);
             if (i > System.Math.PI)
             {
                 i -= (2.0 * System.Math.PI);
             }
             return new Complex(System.Math.Log(r), i);*/
        }

        ///<summary>Given two complex types return the one with the maximum norm</summary>
        public static Complex Max(Complex v1, Complex v2)
        {
            if (Norm(v1) >= Norm(v2))
            {
                return v1;
            }
            else
            {
                return v2;
            }
        }

        ///<summary>Return the Euclidean norm of a complex type</summary>
        public static double Norm(Complex value)
        {
            return ComplexMath.Absolute(value);
        }

        ///<summary>Return the polar representation of a complex type</summary>
        public static Complex Polar(Complex value)
        {
            return new Complex(ComplexMath.Absolute(value), System.Math.Atan2(value.Imaginary, value.Real));
        }

        ///<summary>Returns the sine of a complex type</summary>
        public static Complex Sin(Complex value)
        {
            if (value.IsNaN())
            {
                return Complex.NaN;
            }
            return new Complex(System.Math.Sin(value.Real) * System.Math.Cosh(value.Imaginary), System.Math.Cos(value.Real) * System.Math.Sinh(value.Imaginary));
        }

        ///<summary>Returns the hyperbolic sine of a complex type</summary>
        public static Complex Sinh(Complex value)
        {
            if (value.IsNaN())
            {
                return Complex.NaN;
            }
            if (value.Imaginary == 0)
            {
                return new Complex(System.Math.Sinh(value.Real));
            }
            if (value.Real == 0)
            {
                return new Complex(0, System.Math.Sin(value.Imaginary));
            }
            return new Complex(System.Math.Sinh(value.Real) * System.Math.Cos(value.Imaginary), System.Math.Cosh(value.Real) * System.Math.Sin(value.Imaginary));
        }

        ///<summary>Returns the square root of a complex type</summary>
        public static Complex Sqrt(Complex value)
        {
            Complex result;
            if (value.IsNaN())
            {
                return Complex.NaN;
            }

            if (value.IsInfinity())
            {
                return Complex.Infinity;
            }

            if ((value.Real == 0.0) && (value.Imaginary == 0.0))
            {
                result = Complex.Zero;
            }
            else
            {
                double x, y, t;

                x = System.Math.Abs(value.Real);
                y = System.Math.Abs(value.Imaginary);
                if (x >= y)
                {
                    t = y / x;
                    t = System.Math.Sqrt(x) * System.Math.Sqrt(0.5 * (1.0 + System.Math.Sqrt(1.0 + t * t)));
                }
                else
                {
                    t = x / y;
                    t = System.Math.Sqrt(y) * System.Math.Sqrt(0.5 * (t + System.Math.Sqrt(1.0 + t * t)));
                }

                if (value.Real >= 0.0)
                {
                    result = new Complex(t, value.Imaginary / (2.0 * t));
                }
                else if (value.Imaginary >= 0.0)
                {
                    result = new Complex(value.Imaginary / (2.0 * t), t);
                }
                else
                {
                    result = new Complex(-value.Imaginary / (2.0 * t), -t);
                }
            }

            return result;
        }

        ///<summary>Returns the tangent of a complex type</summary>
        public static Complex Tan(Complex value)
        {
            if (value.IsNaN())
            {
                return Complex.NaN;
            }

            double real = 2.0 * value.Real;
            double imag = 2.0 * value.Imaginary;
            double denom = System.Math.Cos(real) + System.Math.Cosh(imag);

            return new Complex(System.Math.Sin(real) / denom, System.Math.Sinh(imag) / denom);
        }

        ///<summary>Returns the hyperbolic tangent of a complex type</summary>
        public static Complex Tanh(Complex value)
        {
            if (value.IsNaN())
            {
                return Complex.NaN;
            }

            double real = 2.0 * value.Real;
            double imag = 2.0 * value.Imaginary;
            double denom = System.Math.Cosh(real) + System.Math.Cos(imag);

            return new Complex(System.Math.Sinh(real) / denom, System.Math.Sin(imag) / denom);
        }

        /// <summary>Returns the inverse sine of a complex type.</summary>
        public static Complex Asin(Complex value)
        {
            if (value.IsNaN())
            {
                return Complex.NaN;
            }
            Complex ni = new Complex(0, -1);
            return ni * ComplexMath.Log(Complex.I * value + ComplexMath.Sqrt(1 - value * value));

        }

        /// <summary>Returns the inverse cosine of a complex type.</summary>
        public static Complex Acos(Complex value)
        {
            if (value.IsNaN())
            {
                return Complex.NaN;
            }
            Complex ni = new Complex(0, -1);
            return ni * ComplexMath.Log(value + Complex.I * ComplexMath.Sqrt(1 - value * value));
        }

        /// <summary>Returns the inverse tangent of a complex type.</summary>
        public static Complex Atan(Complex value)
        {
            if (value.IsNaN())
            {
                return Complex.NaN;
            }
            return (1 / (Complex.I * 2)) * (ComplexMath.Log((1 + Complex.I * value)) - ComplexMath.Log((1 - Complex.I * value)));
        }

        ///<summary>Returns the inverse hyperbolic sine of a complex type</summary>
        public static Complex Asinh(Complex value)
        {
            if (value.IsNaN())
            {
                return Complex.NaN;
            }
            return -ComplexMath.Log(ComplexMath.Sqrt(1 + value * value) - value);
        }

        ///<summary>Returns the inverse hyperbolic cosine of a complex type</summary>
        public static Complex Acosh(Complex value)
        {
            if (value.IsNaN())
            {
                return Complex.NaN;
            }
            return 2 * ComplexMath.Log((ComplexMath.Sqrt(((value + 1) / 2)) + ComplexMath.Sqrt((value - 1) / 2)));
        }

        ///<summary>Returns the inverse hyperbolic tangent of a complex type</summary>
        public static Complex Atanh(Complex value)
        {
            if (value.IsNaN())
            {
                return Complex.NaN;
            }
            return .5 * (ComplexMath.Log(1 + value) - ComplexMath.Log(1 - value));
        }
    }
}
